function [ Fst, Fr, lamr, ebr ] = factor_st( x, nf, ts, xlags , Flags)
%Calculates structural factors using narrative identification

%Inputs: x - information variables (TxN matrix)
%        nf - number of factors (scalar)
%        ts - identifying dates (vector with dimension <= nf)
%        xlags - number of lags in AR term (scalar)


%Outputs: Fst - Structural factors
%         Fr - Reduced form factors based on PC analysis
%         lamr - Reduced form factor loadings
%         Unexplained variation in the data


%Extract reduced form factors
[Fr,lamr,ebr] = tmp(x,[],nf,0);

%Housekeeping
[~,N] = size(x);
xi = NaN(size(x));
lfx0 = zeros(size(x,1),xlags+1,N);

if length(ts)>nf 
    error('more identifying conditions than factors')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Taking out dynamic component of x's following Stock and Watson (2002) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Setup

    for n = 1:N  
       lfx0(:,1,n) = ones(size(x,1),1);     
            for i = 1:xlags   
                lfx0(:,1+i,n) = lagn(x(:,n),i);   
            end
       
        bb = olss(x(xlags+1:end,n),[lfx0(xlags+1:end,:,n) Fr(xlags+1:end,:)]);
    
         xi(xlags+1:end,n) = x(xlags+1:end,n) - lfx0(xlags+1:end,:,n)*bb(1:xlags+1,:);
    
end
      

%Iteration

if xlags >0
for ii = 1:20;
    
    
      [Fr,~,ebr] = tmp(xi(xlags+1:end,:),[],nf,0);
    
  for n = 1:N;     

      bb = olss(x(xlags+1:end,n),[lfx0(xlags+1:end,:,n) Fr(1:end,:)]);
      xi(xlags+1:end,n) = x(xlags+1:end,n) - lfx0(xlags+1:end,:,n)*bb(1:xlags+1,:);
  end
end
end


%%%%%Factor identification



xip = Fr*lamr';


%Calculate dynamic factor model
lfy0 = ones(size(Fr,1),1);
for i = 1:Flags
    lfy0 = [lfy0 lagn(Fr,i)];
end
fy0v  = Fr(Flags+1:end,:);
lfy0 = lfy0(Flags+1:end,:);

% VAR for DFM
b0   = olss(fy0v,lfy0);   
e0   = fy0v - lfy0*b0;
ex = e0*lamr';

Fst = zeros(size(Fr));

if isempty(ts)==0
    for t = 1:size(Fr,1)
     Fst(t,1:length(ts)) = ex(ts-xlags-Flags,:)'\xip(t,:)';    
    end
end




if length(ts)<nf;

    xi_res = xip-Fst(:,1:length(ts))*ex(ts-xlags-Flags,:);
    
     [F0_res,~,~] = tmp(xi_res,[],nf-length(ts),0);
    
     
     
    Fst(:,length(ts)+1:end)=F0_res;


end

